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^ ■ Abstract 
(N 

^ [ Measuring the similarity of curves is a fundamental problem arising in many application fields. There 

has been considerable interest in several such measures, both in Euclidean space and in more general 
setting such as curves on Riemannian surfaces or curves in the plane minus a set of obstacles. However, 
so far, efficiently computable similarity measures for curves on general surfaces remain elusive. This 
^Cqj [ paper aims at developing a natural curve similarity measure that can be easily extended and computed 

for curves on general orientable 2-manifolds. Specifically, we measure similarity between homotopic 
curves based on how hard it is to deform one curve into the other one continuously, and define this 
"hardness" as the minimum possible surface area swept by a homotopy between the curves. We consider 
I cases where curves are embedded in the plane or on a triangulated orientable surface with genus g, and 

tyj ' we present efficient algorithms (which are either quadratic or near linear time, depending on the setting) 

O ' for both cases. 



> ! 1 Introduction 

l> 
^1 ■ 

Measuring cui^ve similarity is a fundamental problem arising in many application fields, including graph- 
i ics, computer vision, and geographic information systems. Traditionally, much research has been done on 

^ I comparing curves embedded in the Euclidean space. However, in many cases it is natural to study curves 

O ■ embedded in a more general space, such as a terrain or a surface. 

In this paper, we study the problem of measuring curve similarity on surfaces. Specifically, given two 
simple homotopic cui^ves embedded on an orientable 2-manifold (including the plane), we measure their 
similarity by the minimum total area swept when deforming one curve to the other (the "area" of the homo- 
^ . topy between them), and present efficient algorithms to compute this new measure. 



Related work. From the perspective of computational geometry, the most widely studied similarity mea- 
sures for curves is the Frechet distance. Intuitively, imagine that a man and his dog are walking along two 
paths with a leash between them. The Frechet distance between these two paths is the minimum leash 
length necessary for them to move from one end of the paths to the other end without back- tracking. Since 
the Frechet distance takes the "flow" of the curves into account, in many settings it is a better similaiity 
measure for curves than alternatives such as the Hausdorff distance HHH. 

Given two polygonal curves P and Q with n total edges in IR'^, the Frechet distance can be computed 
in 0{n^ log n) time lH. An il{n log n) lower bound for the decision problem in the algebraic computation 
tree model is known ifTTI . and Alt has conjectured that the decision problem is 3SUM-Hard [i2J. Recently, 
Buchin et al. iTTlll show that there is a real algebraic decision tree to solve the Frechet problem with sub- 
quadratic depth, suggesting that perhaps this is not the case. They also give an improved algorithm which 
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runs in 0(n^-v/Iogn(log log n)^) time. Very recently, Agarwal et al. present a novel approach to compute 
the discrete version of the Frechet distance between two polygonal cui"ves in sub-quadratic time |[T1. This is 
the first algorithm for any variant of the Frechet distance to have a sub-quadratic running time for general 
curves. No previous algorithm, exact or approximate, with running time o(n^) is known for general curves, 
although sub-quadratic approximation algorithms for special families of curves are known ll6l 1711231. 

While the Frechet distance is a natural curve similarity measure, it is sensitive to outliers. Variants of 
it, such as the summed-Frechet distance, and the partial Frechet similarity, have been proposed |[T4] [1511241 . 
usually at the cost of further increasing the time complexity. 

The problem of extending and computing the Frechet distance to more general metric space has also 
received much attention. Geodesic distance between points is usually considered when the underlying do- 
main is not M"'. For example, Maheshwari and Yi ll29l computed the geodesic Frechet distance between 
two polygonal paths on a convex polytope in roughly 0{r?K'^\o^{Kv))) time, where n and K are the 
complexity of the input paths and of the convex polytope, respectively. Raichel and Har-Peled consider ap- 
proximating the weak Frechet distance between simplicial complexes in 127 il . Geodesic Frechet distance 
between polygonal curves in the plane within a simple polygon has also been studied jSl [1911251 . 

Rather than comparing distance between only two curves, Buchin et. al. |[T3l propose the concept of a 
median in a group of curves (or trajectories, in their setting). They give two algorithms to compute such 
a median. The first is based simply on the concept of remaining in the middle of the set of curves; this 
algorithm, while fast and simple, has a drawback in that the representative cui^ve might not capture relevant 
features shared by a majority of the input curves. Their second algorithm addresses this issue by instead 
isolating a subset of relevant curves which share the same homotopy type with respect to obstacles that are 
placed in empty regions of the blame; it then computes a medial curve from this relevant subset. 

One issue with generalizing Frechet distance directly to surfaces is that the underlying topology is not 
taken into account; for example, in geodesic Frechet distance, while the length of the leash varies continu- 
ously, the actual leash itself does not. As a result, several measures of similarity have been proposed which 
take the underlying topology into account. Chambers et al. |[T6l proposed the so-called homotopic Frechet 
distance and gave a polynomial (although not efficient) algorithm for when the curves reside in a planar 
domain with a set of polygonal obstacles. The extra requirement for this homotopic Frechet distance is that 
the leash itself and not just its length has to vary in a continuous manner, essentially restricting the homotopy 
class which the leash is in. A stronger variant called isotopic Frechet distance has also been proposed and 
investigated, although no algorithms at all ai^e known to even approximate this distance |[T7l . 

Orthogonal to homotopic Frechet distance is the concept of the height of a homotopy; instead of mini- 
mizing the maximum leash length, this measure views the homotopy as tracing a way for the first curve to 
deform to the second curve, where the goal is to minimize the longest intermediate curve length. Introduced 
independently in two very different contexts |[T0l[T8l . it is not even known if the problem is in NR 

Recent work on approximating the homotopy height and the homotopic Frechet distance has yielded 
efficient 0(log n) approximation algorithms for both of these problems |[26l . However, exact algorithms on 
surfaces for either problem are still unknown. 

New work. In this paper, we develop a natural similarity measure for curves on general surfaces that can 
be computed both quickly and exactly. Intuitively, we measure distances between homotopic curves based 
on how hard it is to deform one curve into the other one, and define this "hardness" as the minimal total 
surface area swept by a homotopy between them, which we call the optimal homotopy area. Our similarity 
measure is natural, and robust against noise (as the area in a sense captures average, instead of maximum, 
deviation from one curve to the other). To the best of our knowledge, this is the first similarity measure for 
curves on general surfaces with efficient polynomial-time algorithms to compute it exactly. 

It is worth noting that this definition in a way combines homotopic Frechet distance with homotopy 
height; those measures compute the "width" and "height" of the homotopy, while our measure calculates 
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the total area. It is thus interesting that while no exact algorithms are known for either of those measures on 
surfaces, we ai^e able to provide a polynomial running time for computing the area of a homotopy. 

We consider both cases where curves are embedded in the plane, or on a closed, triangulated orientable 
surface with genus g. For the former case, our algorithm runs in 0(n log n + P log /) time, where n is the 
total complexity of input curves and / is the number of intersections between them. On a surface, if the input 
is a triangulation of complexity A^, then our algorithm runs in time 0{P log I + ng log n + N). While our 
similarity measure is more expensive to compute for the case of curves in the plane than the Frechet distance 
when / = uj{n), one major advantage is that this measure can be computed on general orientable surfaces 
efficiently. In fact, the ideas and algorithms behind the planar case form the foundation for the handling of 
the case on general surfaces. 

The main ideas behind our approach are developed by examining some properties of one natural class 
of homotopies, including a relation with the winding number of a closed curve. Specifically, the use of the 
winding number enables us to compute the optimal homotopy area efficiently in the plane, where the homo- 
topy is restricted to be piecewise differential and regular. This forms the basis of our dynamic programming 
framework to compute similarity between curves in the plane. We also show how to build efficient data 
structures to keep the total cost of the dynamic program low. 

For the case where the underlying surface is a topological sphere, we extend the winding number in a 
natural way and show how to adapt our planai^ algorithm without additional blow-up in the time complexity. 
For the case when the surface has non-zero genus, we must extend our algorithm to run efficiently in the 
universal cover (which is homeomorphic to the plane) by using only a small portion of it. 

We remark that the idea of measuring deformation areas has been used before in practice Il20ll30l . For 
example, similarity between two convex polygons can be measured by their symmetric difference |l3] |36] ; 
we note that this is not equivalent to homotopy area, although it may be the same value in some situations. 
In another paper, the area sandwiched between an x-monotone curve and another cui-ve is used to measure 
their similarity Q. However, computing the "area" between general curves has not been investigated prior 
to this work. 

2 Definitions and Background 

Paths and cycle. We will assume that we are working on an orientable 2-manifold M (which could be the 
plane). A curve (or a path) on a surface M is a map P : [0, 1] — )• M; a cycle (or a loop) is a continuous 
map 7 : 5^ — )■ Af where is the unit circle. A curve P or a cycle 7 is simple if P{ti) 7^ P{t2) (resp. 
7(ti) / 7(t2)) for any ti / t2. 

Homotopy A homotopy between two paths P and Q (with the same endpoints) is a continuous map H : 



[0, 1] X [0, 1] ^ M where H{0, •) = P, H{1, •) = Q, H{-,0) = P(0) = Q(0) and H{-, 1) = P(l) = 



A homotopy describes a continuous deformation between the two paths or curves: for any value t E [0, 1], 
we let Ht = H{t, ■) be the intermediate curve at time t, where Hq = P and Hi = Q. 

We define the area of a homotopy H to be the total area covered by the image of the homotopy on the 
surface, where an area that is covered multiple times will be counted with multiplicity. More precisely, given 
a homotopy H whose image is piecewise differentiable. 



. The minimum homotopy area between P and Q is the infimum of the areas of all homotopies between P 
and Q, denoted by a {P,Q). If such an infimum does exist and can be achieved by a homotopy, we call that 
homotopy an optimal homotopy. 
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We note that it is not immediately clear that this value exists, depending on the cui^ves and underlying 
homotopy. Minimum area homotopies were considered by Douglas |[22l and Rado [i32i in the context of 
Plateau's problem; they noted that not only is the integral improper in general, but the infimum itself may 
not be continuous. The eventual proof that these exist in reUes on a definition using Dirichlet integrals 
which ensure (almost) conformal parameterizations of the homotopy. See the book by Lawson ll28l for an 
overview of this result as well as several extensions to minimal area submanifolds in more general settings. 

However, beyond a proof of existence, we are interested in computing such 
homotopies, or at least measuring their actual area, in much simpler settings such 
as or a surface. To this end, we restrict the input curves to be simple curves 
which consist of a finite number of piecewise analytic components. We also need H to be continuous and 
piecewise differentiable, so that the integral can be defined. Finally, we will also require that at any time 
t, the intermediate cui^ve Ht is regular (see ll37l for smooth curves and IIBTI for piecewise-linear curves). 
Intuitively, this means that the deformation is "kink"-free OTI . and cannot create or destroy a local loop as 
shown in the right figure (the singular point in the right curve is a kink). Note that this is required for the 
minimum homotopy to even exist; again we refer the reader to the book by Lawson ll28l for details. 

Decomposing arrangements. Consider two simple piecewise analytic curves P and Q with the same 
endpoints. Their concatenation forms a (not necessarily simple) closed curve denoted by C = P o rev(Q), 
where rev{Q) is the reversal of Q. Let Arr(C) denote the arrangement formed by C, where vertices in 
Arr(C) are the intersection points between P and Q. An edge / arc in Arr(C) is a subcurve of either P or 
Q. 

We give C (and thus P and Q) an arbitrary orientation. Hence we can talk 
about the sidedness with respect to C at a point p ^ P. Specifically, a point x G 
is to the right of C at p if it is a counter-clockwise turn from the orientation of the 
vector px the orientation of (tangent of) C at p (see the right figure for an example). 
Given two oriented curves 71 and 72, an intersection point p of them is positive if it 
is a counter-clockwise turn from the orientation of 71 to that of 72 at p. For a curve 7 and a point x G 7, the 
index ofx is the parameter of x under the arc-length parameterization of 7. We sometimes use x to represent 
its index along 7 when its meaning is clear from the context. Given two points j;, y G 7, we will use 7[x, y] 
to denote the unique sub-curve of 7 between points x and y. 

We say that a homotopy H from P to Q is right sense-preserving if for any t,s £ [0, 1], we have that 
either Ht^fit{s) = Ht{s) or Ht+dt{s) is to the right of the oriented curve Ht at Ht{s). If it is the former 
case, then we say that p = Ht{s) is & fixed point at time t. Similarly, we say that H is left sense-preserving 
if for any t, s G [0, 1], Ht{s) is either a fixed point or deforms to the left of the curve Ht. Our homotopy 
H is sense-preserving if it is either right or left sense-preserving. The sense-preserving property means that 
we can continuously deform the curve P always in the same direction, without causing local folds in the 
regions swept. Intuitively, any optimal homotopy should have this property to some extent, which we will 
make more precise and prove later. 

3 Structure of Optimal Homotopies 

Given two simple curves P and Q (with the same end points) embedded on an orientable 2-manifold M, 
let X = {xi, . . . , X/} denote the set of I intersection points between them, sorted by their order along P. 
Given a homotopy H from P to Q, a point p G Af is called an anchor point with respect to H if it remains 
on H{t, ■) = Ht at all times t G [0, 1]. Of course not all intersection points are anchor points. However, 
if p is an anchor point, then it is necessarily an intersection point between P and Q, as p G Hq = P 
and p £ Hi = Q. We exclude the beginning and ending end points of P and Q from the list of anchor 
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points, as they remain fixed for all homotopies. In what follows, we show that any optimal homotopy can 
be decomposed by anchor points such that each of the resulting smaller homotopies has a simple structure. 

Specifically, consider an arbitrary optimal homotopy H*. Let B = {bi, . . . , b^} be the set of anchor 
points with respect to H*, the minimum area homotopy. We order the bj's by their indices along P. It 
turns out that the order of their indices along Q is the same, and the proof of this simple observation is in 
Appendix lAl 

Observation 3.1 The order of hi 's along P and along Q are the same. 

This observation implies that we can decompose H* into a list of sub-homotopies, where H* moiphs 
P[bj,bj+i] to Q[bj,bj+i]. Obviously, each H* is necessarily optimal, and it induces no anchor points. 
The following result states that an optimal homotopy without anchor points has a simple structure, which is 
sense-preserving. Intuitively, if any point changes its deformation direction at any moment, the defomiation 
will sweep across some area redundantly and thus cannot be optimal. The detailed proof is in Appendix |B] 

Lemma 3.2 If an optimal homotopy H from P to Q has no anchor points, then it is sense-preserving. 

4 Minimum Area Homotopies In The Plane 

In this section, we consider the case where the input consists of two simple polygonal curves in the plane. 
We develop an algorithm to compute the similarity between P and Q in 0{I^ log I + n log n) time, where 
n is the total complexity of input curves and / is the number of intersections. Note that / = Q{n?) in 
the worst case, although of course it may be much smaller in some cases. Although efficient algorithms 
for comparing curves in the plane exist (such as the Frechet distance), our planar algorithm will be the 
fundamental component for comparing curves on general surfaces in the next section. It turns out that our 
approach can easily be extended to measure similarity between simple cycles in the plane; see Appendix IdI 

4.1 Relations to Winding Numbers 

We are given two simple (open) curves in the plane which share common endpoints. Previously, we have 
shown that if an optimal homotopy does not induce anchor points, then it is sense-preserving. The impli- 
cation of this result is manifested by using the winding number, defined for a loop 7 in the plane at a base 
point X. 

Intuitively, imagine starting from a point y on 7, and connecting x and y by a string. The winding 
number at x w.r.t. 7, denoted by wn(x; 7), is an integer measuring how many times this string winds, in a 
clockwise manner, around x as y traverses 7. More formally, consider an infinite ray / based at x which is 
generic (so it has a finite set of transversal intersections / crossings with 7). Consider a crossing 7(t) between 
the ray / and 7. This crossing is positive if the triangle x, 7(t), and 7(t + e) is oriented counterclockwise, 
and is negative if oriented clockwise. Then wn(a;; 7) is the number of positive crossings minus the number 
of negative crossings with respect to any generic ray from x. 

We say an oriented curve 7 has consistent winding numbers if 
wn(a;, 7) is either all non-negative, or all non-positive, for all x G IR^. 
Note that for a curve with consistent winding numbers, we can always 
orient the curve appropriately so that wn(x, 7) is all non-negative. Two 
examples are shown in the figure on the right, where the second exam- 
ple has consistent winding numbers. Let Arr(7) denote the arrangement formed by the curve 7. All points 
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in the same cell of the aixangement Arr(7) of 7 have the same winding number, and the winding num- 
bers of two neighboring cells differ by 1. The relation of consistent winding numbers and sense-preserving 
homotopies is given below, and the proof can be found in Appendix ??. 

Lemma 4.1 If there is a sense-preserving homotopy H from P to Q, then the closed curve P o rev{Q) has 
consistent winding numbers. 

Proof: Without loss of generality, assume that the map H is left sense-preserving, always deforming an 
intermediate curve to its left. Consider the time-varying function F : [0, 1] x IR^ — )• Z, where F{t, x) = 
wn(a;; Ht) is the winding number at x G with respect to the curve parameterized by Hf. Obviously, 
F(0, x) = wn(x; P o iev{Q)), and x) = 0. During the deformation, x) changes by either 1 or 
— 1 whenever the intermediate curve sweep over it. Since the homotopy is left sense-preserving, when an 
intermediate curve sweeps x, x always moves from the left side of the intermediate curve to its right side. 
Hence the winding number x decreases monotonically. Since in the end, the winding number at each point 
is zero, wn(x; P o rev((5)) = F{0, x) > 0. 

If the map H is left sense-presei^ving, then a symmetric argument shows that wn(2;; Porev((5)) < for 
all X e ]R2. ■ 

Next, we describe two results to connect the above lemma to the computation of optimal homotopy. 
First, we define the total winding number Tw(7) of a curve 7 as 

Tw(7) = / wnfx; 7)(iz^(x), 

where di'{x) is the area formQ. The following observation is straightforward. 
Observation 4.2 For any P and Q in the plane, 

a{P,Q) > |Tw(Porev(Q))|. 



Proof: Take any regular homotopy H from P to Q. The area of a regular homotopy H in our setting can be 
reformulated as an integral on the image domain as 

Area{H) = / degH{x)di'{x), 

where degnix) is defined as the number of connected components in the pre-image of x under H. In other 
words, degnix) is the number of times that any intermediate curve Ht sweeps thi^ough x. Now consider the 
function F : [0, 1] ^ IR defined as F{t) = Tw{H^ o Q). Obviously, F(0) = Tw(P o Q), F(l) = 0, and 
F is a continuous function. Furthermore, each time the winding number at a point x changes by 1 for some 
t G [0, 1], it means that some intennediate cui've H{t) sweeps thi^ough it. Hence |wn(x)| is a lower bound 
for deguix). We thus have that 

|Tw(Porev(Q))| < / \wTi{x)\dv{x) < / degH{x)dv{x) 

for any regular homotopy H, implying that 

|Tw(PoQ)| < \niArea{H) = a{P,Q). 



'Note that this allows us to use any Riemannian metric on the plane (including the standard Euclidean metric). This will be 
necessary later when we use the same algorithm for curves in a universal covering space. 
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(a) (b) (c) (d) 

Figure 1 : (a) The cell R with highest positive winding number. Its boundary consists of alternating P-arcs (red) and 
Q-arcs (green). The two cases of relations between P[p, q] and R are shown in (b) and (d), respectively. For case (b), 
we can deform P to sweep through il as shown in (c), and reduce the number of intersections by 2. 

Lemma 4.3 Given P and Q, if j = P orev[Q) has consistent winding numbers, then a{P, Q) = |Tw(7)|. 



Proof: We provide a sketch of the proof here to illustrate the main idea; see [?] for the full proof. We prove 
the claim by induction on the number of intersections between P and Q. The base case is when there is no 
intersection between P and Q. In this case, 7 is a Jordan curve which decomposes into two regions, one 
inside 7 and one outside. The optimal homotopy area a{P, Q) in this case is the area of the bounded cell. 
All points in the bounded cell have winding number 1 (or —1) and the claim follows. 

Now assume that the claim holds for cases with at most k—1 intersections. We next prove it for the case 
with k intersections. Let an X-arc denote a subcurve of curve X. Consider the arrangement Arr(7) formed 
by 7 = Porev((5). Since P and Q are simple, every cell in this arrangement has boundary edges alternating 
between P-arcs and Q-arcs. Assume without loss of generality that 7 has all non-negative winding numbers. 
Consider a cell R e Arr(7) with largest (and thus positive) winding number. Since its winding number is 
greater than all its neighbors, it is necessary that all boundary arcs are oriented consistently as shown in 
Figure [T](a), where the cell R (shaded region) lies to the left of its boundary arcs. 

If R has only two boundary arcs, e from P and e' from Q, respectively, then we 
can moiph P to another simple curve P' by deforming e through R to rev{e') as 
illustrated on the right. The area swept by this deformation is exactly the area of cell 
R. Furthermore, after the deformation, every point x ^ R decreases their winding 
number by 1, and no other point changes its winding number. Since any point in this 
cell initially has strictly positive winding number, the resulting curve -)' = P' o Q 
still has all non-negative winding numbers. The number of intersections between P' and Q is /c — 2. By 
induction hypothesis, a{P',Q) = Tw(7'). Since Tw(7) - Tw(7') = Area(R), we have that Tw(7) = 
a{P', Q) + Area(R). It then follows from Observation O and the fact a{P, Q) < a{P', Q) + Area(R) 
thata(P,Q) = Tw(7). 

Otherwise, the cell R has more than one P-arc. Take the P-arc ei with the smallest index along P, 
and let p be its ending point. Let 62 be the next P-arc along the boundary of R, and q its starting point, 
and Q[q,p] the Q-arc between ei and 62, denoted by e in Figure [T] P\p,q] and —Q\p,q\ bound a simple 
polygon, which we denote by Q.. Since Q. does not intersect R, either Q. is on the opposite side of the Q-arc 
e from the inteiior of R (Figure [T](b)), or they are on the same side (Figure [T](d)). It turns out that in both 
cases, we can deform P to another simple curve P' such that (i) the number of intersections is reduced by 
2, and (ii) P' o Q still has consistent winding numbers. For example, in the case of Figure [r|(b), P is then 
deformed to sweep the region Q, as shown in Figure [T](c). By applying the induction hypothesis to P' o Q, 
we are able to obtain that a{P,Q) = Tw(7). The details are in Appendix O ■ 




7 



4.2 The Algorithm 



Lemma |3^ and |4~T] imply that if the closed curve P o rev((5) produces both positive and negative winding 
numbers, then any optimal homotopy from P to Q must have at least one anchor point. On the other hand, if 
it has consistent winding numbers, then by Lemma 1431 we can compute the optimal cost to deform them by 
simply computing the total winding number. This leads to a simple dynamic-programming (DP) approach 
to compute a{P, Q). 

Specifically, let xq, xi, . . . , x/ denote the intersection points between P and Q, ordered by their indices 
along P, where xq and x/ are the beginning and ending points of P and Q, respectively. Let T{i) be 
the cost of the optimal homotopy between P[0,Xj] and Q[0,Xj], and C[i,j] the closed curve formed by 
P[xj, Xj] o Q[xj, Xj]. We say that a pair of indices (z, j) is valid if (1) x, and Xj have the same order along 
P and along Q; and (2) the closed curve C[i,j] has consistent winding numbers. We have the following 
recursion: 



{Tw(C[0, z]), If C[0, i] has consistent winding numbers; 

i^^<., (j, i) is valid i Tw(C[i, i]) + T{j) }, Otherwise. 



4.3 Time Complexity Analysis 

The main components of the DP framework described above are (i) to compute Tw(C[i, j]) for all pairs of 
js, and (ii) to check whether each pair (i, j) is valid or not. These can be done in 0{I^n) total time in a 
straightforward manner. We now show how to compute them in 0(/^ log /) time after 0{I log I + n log n) 
pre-processing time. Specifically, we describe how to compute such information in 0(/log/) time for all 
C[r, i]s for a fixed r G [1, /] and all indices i > r. 

To simplify the description of the algorithm, we extend Q on both sides until infinity, and obtain Q. Now 
collect all intersection points between P and Q, {xi, . . . , x/}, which is a super-set of previous intersection 
points, and sort them by their order along the curve P. The algorithm can be made to work with Q directly, 
but using Q makes the intuition behind our algorithm, as well as its description, much more clear. 

Note that Q divides the plane into two half-planes. For illustration purpose, we will draw Q as a hor- 
izontal line, and use the upper and lower half-planes to refer to these two sides of Q. Another way to see 
that regarding Q as a horizontal line does not cause any loss of generality is that one can always find a 
homeomorphism from — IR^ such that the image of Q is a horizontal line under this homeomorphism. 

Now for a fixed integer r, we traverse P starting from x^. We aim to maintain appropriate data structures 
so that each time we pass through an intersection point Xj with Q, we can, in 0(log/) time, (1) check 
whether (r, i) is valid, and (2) obtain total winding number for C[r, i]. 




Figure 2: Illustration of the regions R^s. 



Total winding numbers. We first explain how to maintain the total winding number for the closed curve 
C[r, i] = P[xr.,Xj] o (5[xj,Xr] as i increases. Assume i changes from u to u + I. Since x„ and x„-|_i are 
two consecutive intersection points along P, the arcs P[x„,x„+i] and Q[^u,^u+i] form a simple closed 
polygon which we denote by i?„ (shaded region Figure |2]i. Comparing the aiTangement Avr {C[r,u + 1]) 
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with Arr (C[r, u]), regardless of where r is, only points within will change their winding number, either 
all by +1 or all by —1, depending on whether i?„ is to the right side or the left side of the P-arc P[^u, ^u+i], 
respectively. The winding numbers for points outside Ru are not affected. Hence the change in the total 
winding number is simply Q„Area(Ru), where q„ is either +1 or —1. See Figure|2j where all points in Ru 
will decrease their winding number by 1 as we move from C[r, u] to C[r, u + 1]. 

We can pre-compute the area of R^s for all n in 0(n log n + / log /) time, by obsei^ving that the set of 
RuS satisfy the parentheses property: Namely, either R^ and R^ are, disjoint in their interior, or one contains 
the other. 



p 




Figure 3: The containment relations of all i?„ regions can be represented as a forest structure on the right. 

Specifically, first, we compute the arrangement of Aii{P+Q) and the area of all cells in it in 0{n log n+ 
/) time. Each Ru is the region bounded between a P-arc P[xj,Xj+i] and a corresponding Q-segment 
Q[xj,Xj+i]. Since no two P-arcs intersect, the containment relationship between such P-arcs satisfies 
parentheses property. In particular, we can use a collection of trees to represent the containment relation 
among all regions P„s. See Figure |3]for an illustration. The difference between the region represented at a 
parent node and the union of regions represented by all its children is a cell in Arr(P + Q). For example, the 
shaded cell in the right figure is the difference between Rq and its children R2 and P4. We can thus compute 
the area of all P„s by a bottom-up traversal of these trees. Computing these trees take 0(/log/) time by 
first sorting all intersection points with respect to their order along Q. Traversing these trees to compute all 
i?„s takes 0(1) time. Putting everything together, we need 0(n log n + / log /) time. 

With the area of all P^s known, updating the total winding number from C[r, u] to C[r,u + 1] takes 
only constant time. 

Checking the validity of (r, i)s. To check whether (r, i) is valid or not, we need to check whether all 
cells in the arrangement Arr(C[r, i]) have consistent winding numbers. First observe that for any r and 
i, Arr(P + Q) is a refinement of the arrangement Arr(C[r, i]). That is, a cell in Arr(P + Q) is always 
contained within some cell in Arr(C[r, z]). Hence all points within the same cell of Ait(P + Q) always 
have the same winding number with respect to any C[r, i]), and we simply need one point from each cell in 
Arr(P + Q) to maintain the winding number for all cells in Arr (C[r, i]), for any r and i. We now describe 
how to maintain the winding number for cells of Att{P + Q) (thus for Arr(C[r, n])s) as we pass each u > r, 
so that we can check whether C[r, u] has consistent winding numbers or not efficiently. 

To this end, take four points around each intersection point Xj of P and Q 
(shown as stars in the right figure). The collection of such representative points 
hit all cells in Arr(P + Q). (It does not matter whether there may be more than 
one point taken from a cell of Arr(P + Q).) Hence Arr(C[r, i\) has consistent 
winding number if and only if all these representative points have consistent 
winding numbers. Next, we build a data structure to maintain the winding numbers for these points as i 
increases. Specifically, let U be the set of representatives that are to the right of Q, which ai^e the stars above 
Q in the right figure. (Those to the left of Q will be handled in a symmetric manner). Each point has a key 
associated with it which is its index along Q. We build a standard balanced 1-D range tree on U based on 
such keys, where each leaf / stores a point from U. Every internal node v is associated with an interval 
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[Z^, r^,], where and r„ are the smallest and lai^gest keys stored in the subtree rooted at v. In other words, 
all representatives with an index along Q within r^] are stored in the subtree rooted at v. At every node 
V, interior or not, we also store a value addWy. To compute the winding number for the representative point 
Pf stored at a leaf node /, we identify the path {^0,^1, ■ ■ ■ ,Va = /} from the root vq to /. The winding 
number for pj is simply Yli=o O'ddWy. . Finally, each internal node v also stores the maximum and minimum 
winding numbers associated with all leaves in its subtree. At the beginning, all winding numbers ai^e zero. 
The size of this tree is 0{I) with height 0(log /), and can be built in 0{I log I) time once the arrangement 
Arr(P + Q) is known. 

Let qi denote the index of point Xj along Q (or can be considered as the x-coordinate of Xj). As we move 
from C[r, u] to C[r, u + 1], cells of Arr(P + Q) contained in Ru should either all increase or all decrease 
their winding number by 1. Note that representatives of these cells are simply those contained in the interval 
[Qm, qu+i] (or [q^+i, qu] if q^+i < q^). Hence updating the winding number is similar to an interval query 
of [q„, qu+i], and the 0(Iog /) number of nodes in the canonical decomposition of [qu, qu+i] update their 
addW^ values by either +1 or —1 depending on the sideness of Ru with respect to the arc P[x„,x„+i]. 
The minimum and maximum winding numbers can also be updated 0(1) time per visited node. The entire 
process visits O(logl) nodes, and thus takes 0(log /) time as i increases from u to n + 1. To see whether 
C[r,u + 1] has consistent winding numbers or not, we only need to check the minimum and maximum 
winding numbers stored at the root of the tree, denoted by Wmin and Wmax, respectively. If Wmin x Wmax 
equals to zero, then all winding numbers w.r.t. C[r,u + 1] are either all non-negative or all non-positive. 
Otherwise, (r, u + 1) is not valid. 

Repeat the above process for every r G [1, Overall, after 0{{n + I) log n) pre-processing, we can 
check whether (r, i) is valid or not and compute Tw(C[r, i]) for all r S [1, /] and alH > r in 0(/^ log /) 
time. Putting everything together, we have the following result. 

Theorem 4.4 Given two simple polygonal chains P and Q (with the same endpoints) in the plane ofn total 
complexity, and with I intersection points between them, we can compute the optimal homotopy and its area 
in 0{P log I + n log n) time and 0{I^ + n) space. 

The case where we have two simple cycles P and Q in the plane is discussed in Appendix and we 
obtain the following extension: 

Corollary 4.5 Given two polygonal cycles P and Q in the plane ofn total complexity and with I intersec- 
tion points, we can compute the optimal homotopy and its area in 0{I{I^ log / + nlogn)) time if I > 0; 
and compute the optimal homotopy area in O (nlogn) time if I = 0. 

5 Minimum Area Homotopies on 2-ManifoIds 

In this section, we consider optimal homotopy between curves P and Q on an orientable and triangulated 
2-manifold M without boundary. Our input is a triangulation K of M with complexity N, and two simple 
homotopic polygonal curves P and Q sharing endpoints. Edges in P and Q are necessarily edges from the 
triangulation K. The total complexity of P and Q is n, and there are / number of intersections between 
them. Note that in this setting, / = 0(n). Below we discuss sepai'ately the cases when M has non-zero 
genus and when M is a topological sphere. 

5.1 Surfaces with Positive Genus 

Given an orientable 2-manifold M, let U{M) be a universal covering space of M with (f) : U{M) — )• M 
the associated covering map. Note that (f) is continuous, surjective, and a local homeomorphism. (For full 
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details on covering spaces, we refer the reader to topology textbooks that address this area |]33l; we will 
also build on existing algorithmic techniques developed for the computing and working in the universal 
cover lEB [Ml-) 

For any path 7 in M, if we fix the lift (pre-image) of its starting point, then it lifts to a unique path 7 
in U{M), such that (/)(7) = 7. Since P and Q are homotopic with common endpoints, the closed curve 
formed by C = P o Q is contractible on M, and the lift of C, denoted by C, is a closed cui"ve in U{M). By 
the Homotopy Lifting Property of the universal cover Il33l, we have: 

Observation 5.1 Once we fix the lift of the starting point of P and Q in U{M), there is a one-to-one 
correspondence between homotopies between P and Q in M and those between P and Q in U{M). 

We now impose an area measure in U{M) by lifting the area measure in M; this can be done via the 
map (f), which is a local homeomorphism. Now the area of a homotopy in M is the same as the area of its 
lift in U{M). As such, we can convert the problem of finding an optimal homotopy in M to finding one 
in U{M). Furthermore, for any orientable compact 2-manifold with genus g > its universal cover is 
topologically equivalent to M^. Intuitively, this means that we can then apply algorithms and results from 
previous section to the universal covering space. 

More specifically, our algorithm proceeds as follows: 

Step 1: Compute relevant portion of U{M). We will construct a portion of a universal covering space 
U{M) made from polygonal schema of M |[35ll2TI . Specifically, we use the algorithm from 11211 to construct 
a reduced polygonal schema T in 0{N) time. The universal covering space consists of an infinite number 
of copies of this polygonal schema glued together appropriately. We call each copy of the polygonal schema 
in the constructed universal covering space a tile. 

Recall that the universal covering space U{M) is homeomorphic to M^. We fix a lift of the starting 
endpoint of P and Q in U{M) and obtain a specific lift P and Q for P and Q respectively. Since P and 
Q are homotopic, P and Q fomi a closed curve, denoted hy C = P o rev((5). Note that the number of 
intersection points between P and Q is at most /, as every intersection point in the lift necessarily maps to 
an intersection point of P and Q under 0, but not vice versa. 

Consider the arrangement formed by Arr(C') in the planar domain U (M). We will construct the portion 
of the universal covering space U C U{M) which is the union of tiles that intersect or are contained inside 
of Arr(C'). 

From ll2n . we know that the Ufted curve C passes through 0(n) tiles in U{K). However, while the 
total number of tiles in the interior of Arr(C') is 0{n) for the case where 51 > 1, it can be Q{n?) for the 
case when g = \. Hence we will separate the case iox g = 1 and g > 1, since we wish to avoid the O(n^) 
overhead in the genus 1 case. 

For the case (7 > 1, we use the algorithm by Dey and Schipper II2TI to compute the relevant portion 
U of the universal covering space in 0{n\ogg + N) time. The output contains all 0{n) copies of the 
polygonal schema in U, where each tile is represented by a reduced Ag-gon without being explicitly filled 
with triangles from K. However, recall that P and Q are cui^ves which follow edges of the triangulation; 
in this construction of the polygonal schema tiles, each edge of K can be broken into 0{g) pieces. So in 
the worst case, we must break each edge in P or Q into 0{g) pieces, giving a total complexity for P and 
Q is 0{ng) in this representation of U. Once these are known, we can compute the combinatorial structure 
of the arrangement of C in [/, as well as the description of the set of tiles each cell in Arr(C') intersects or 
contains, in 0{ng + I log /) time. 

For the case g = 1, the input manifold is a torus, and the canonical polygon schema for it is a rectangle 
with oriented boundary arcs aha~^h~^. Imagine now that we give the base polygonal schema Tq (which 
is the tile that contains the lift of the starting point of P and Q) a coordinate (0, 0); we must now assign a 
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(a) (b) 

Figure 4: (a) A combinatorial view of the universal covering space U{M). a and h are the generators and we can 
give each cell a coordinate, (b) The lift of P (solid curve) and the lift of Q (dashed curve). The heavily shaded region 
are copies of polygonal schema contained inside cells of Arr(C), and their total number can be easily computed by a 
scanning algorithm. i?i is an essential cell; R2 and are two non-essential cells. 

coordinate for every other copy of the polygon (as shown in Figure Ha)). Specifically, a copy of polygonal 
schema T has coordinate (z, j) if the closed loops whose lifts start in Tq and end in T have the same 
homotopy type as a^V . We can obtain the sequence of the rectangles (and their coordinates) that the curve 
C will pass through in 0(n + N) time Once these coordinates are known, the combinatorial structure 
of the arrangement of (7 in [/ can also be computed in 0(n + / log /) time. Note that in this case, we have 
avoided explicitly enumerating the set of tiles fully enclosed within Arr(C') (the shaded tiles in Figure |4] 
(b)), whose number can be O(n^) instead of 0{n) when g = 1. 

Step 2: Area of cells in Arr(C'). In order to perform our algorithm introduced in Section 5]to the lifted 
curves P and Q, in addition to the combinatorial structure of Arr(C'), we also need the area of each cell in 
Arr(C'). We first describe how to compute it for the case g = I. 

Take any cell X in Arr(C') and assume the boundary of X intersects m copies of polygonal schema. 
Even though that X may contain 0(m^) copies of (rectangular) tiles in its interior, we do not need to 
enumerate these interior tiles explicitly to compute their total area. 

Indeed, by a scanning algorithm from left to right, we can compute in 0{m) time how many tiles are 
completely contained inside X (heavily-shaded regions in Figure |4](b)) (note that the coordinates of each 
tile traversed by the boundary of X are known). Since the area of every polygonal schema is simply the 
total area of the input triangulation, we can compute the total area of tiles contained inside X in 0{m) time. 

Now let R be the collection of tiles that intersect the boundary of X. It remains to compute the total 
area of R n A". Call each region in Tn A a sub-cell, for any tile T G R. Let G denote the boundary curves 
of the polygonal schema T. There are two types of sub-cells: the essential ones which contain at least one 
intersection point between P and Q as their vertices, and the non-essential ones which have no intersection; 
see Figure |4] (b) for examples. Note that a non-essential cell is bounded by arcs from G alternating with 
P/Q-aics from P or Q, since there is no intersection of P and Q along the boundary of a non-essential cell. 
(Here, a P/Q-aic refers to either a P-arc or a Q-aic). 

First let us consider the collection of non-essential sub-cells formed by alternating G-arcs (boundary arcs 
of a tile) and arcs from P and Q, and compute the area of each such non-essential sub-cells. If we plot all 
the P-arcs within a single tile T, no two P-arcs can intersect in this tile, since P is a simple curve. Imagine 
that we pick an arbitrary but fixed point on the boundary G of the polygonal schema T as the origin o. Each 
P-arc a subdivides T into two regions; we let denotes the canonical one excluding o. Note that since 
P is a simple curve, the set of canonical regions T^s for all P-arcs must satisfy the parenthesis property, 
and these regions and their areas, called canonical areas, can be computed in 0{ng log n + N) time using a 
data structure similar to one used in Section l43] to compute the area of RuS. See Figure |5]for an illustration. 
Similar, we can put all Q-arcs within the same tile and compute the canonical regions / areas associated with 
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(a) 



(b) 



Figure 5: (a) We overlay all non-essential sub-cells involving P-arcs into one copy of the polygonal schema, 
(b) An example of the canonical region Tq, is shown for arc a (shaded region in the top-right corner). The 
shaded region in the middle is a sub-cell X which can be computed as T^^ — T^^ — T/jj — T^g, where /JjS 
are the boundary P- and Q-arcs for X. Among these P/Q-aics, /3o is the top-most arc in the containment 
relation. 

all Q-arcs in 0{ng log n + A^) time. Once these areas are known, the area of each non-essential sub-cell 
can be computed in 0{s) time where s is the number of P-arcs and Q-arcs on the boundary of this sub-cell: 
Specifically, it is the difference between the canonical area of the top-most P/Q-aic and the union of the 
canonical areas of all other P/Q-arcs on the boundary of this cell. See Figure |5](b). Hence the areas of all 
non-essential sub-cells can be computed in 0{n) time once all T^s are known. The total time complexity 
required here is thus 0{ng log n + A^). 

What remains is to compute the area of all essential sub-cells. Note that there are 0{I) essential sub- 
cells since each contains an intersection between P and Q. Let a PQ-aic to refer to an ai^c that starts and 
ends with points on G (the boundary of the polygonal schema T) and consists of alternating P- and Q-arcs. 
An essential sub-cell is either completely contained within a polygonal schema, or its boundary consists 
of PQ-arcs, G-arcs, P-arcs and Q-arcs where no two such arcs can be consecutive: they are separated by 
G-arcs. Now collect all P-arcs and Q-arcs that are involved in the boundary arcs of those essential sub-cells 
completely contained within a tile. Plot them within the same tile T and compute their aixangement A as 
well as the area for each cell in A. This can be done in 0{ng log n + A^) time. Since A can have only 0{I) 
vertices in the interior of the tile T, A contains 0{I) cells. If an essential sub-cell X is completely contained 
within a polygonal schema, then it is a union of a set of cells from A. We can simply spend 0{I) time to 
go through cells in A, identify those contained in X and return their total ai^ea. Hence it takes 0{P) time to 
compute the area of all such 0{I) essential sub-cells. If an essential sub-cell X has G-arcs on its boundary, 
then we need a slightly more complicated way to handle it. 

Specifically, for all the remaining essential sub-cells, there can be 0{I) number of PQ-arcs along their 
boundaries, denoted by L. We collect all P-arcs and Q-arcs involved in L and plot them in the same tile 
T and compute their amngement Arr(L). Each PQ-arc a G L divides the tile T into two regions, and 
we define Tq to be the canonical one that excluding a specific origin o on G similar to before. consists 
of a union of cells from the arrangement Arr(L), and we can compute the area of Tq, in 0{I) time since 
Arr(L) has 0{I) cells. Overall, in 0(/^) time, we can compute the area of all T^s for all PQ-ai^cs a G L 
Now take an essential sub-cell X that has s number of P-, Q-, or PQ-arcs along its boundaiy, denoted by 
Q!i, a2, . . . , a^. Let ai be the arc (which can be P-, Q- or PQ-arc) whose endpoints along G spans the 
largest interval. Then, X can be represented as AT = T^^ — IJt6[2 s] -^^i' where Tq,- is the canonical region 
defined by an arc Oj. Since the area of all canonical regions are known (for P-aixs or Q-arcs, we have 
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computed their canonical areas before), X's are a can be computed in 0{s) time. Computing the area of all 
remaining essential sub-cells thus takes 0{P + n) time. 

Putting everything together, the total time needed to compute the area of all cells in Arr((7) is 0{ng log n+ 
N + P) when g = I. The case when g' > 1 is similar but much simpler. Indeed, we now can afford to 
compute all the tiles contained within any cell of Arr(C') explicitly, as their total number is bounded by 
0(?i) lliniMl . The areas of essential and non-essential sub-cells are computed using the same algorithm as 
above. The total time complexity is 0{ng log n + N + l'^). 

Step 3: Putting everything together. With the combinatorial structure of Arr(C') and the area of each cell 
computed, we now apply the algorithm from Section |4]2]to compute the optimal homotopy in 0{I^ log / + 
n^rlogn) time mU{M), which, by 15.11 gives the optimal homotopy between P and Q in M in the same 
time bound. The total time complexity for the entire algorithm is 0{ng log n + log / + N). 

5.2 The Case of the Sphere 

We now consider the remaining case where the input has g = 1, so M is a (topological) sphere S. All paths 
on S are homotopic. The universal cover of a sphere is itself, and hence is compact. However, the previous 
algorithm in Section I42] works for a domain homeomorphic to and cannot be directly applied. We now 
sketch how we handle the sphere case. Missing details can be found in Appendix El 




Figure 6: Two ways of sweeping a curve on sphere from base point p. 

We observe that the results in Section |3] still hold. However, as the sphere is compact, the winding num- 
ber is not well-defined. For example, see Figure |6] where there are two ways that the curve 7 winds around 
the point p. In the first case, the winding number at p is 0, while in the second case, the winding number 
is —1. In order to use a dynamic programming framework as before to compute the optimal homotopy 
between P and Q, we need to develop analogs of Lemma 14.11 and 14. 3 1 for curves on the sphere. 

To this end, note that if we remove one point, say z G S from the sphere S, then the resulting space 
Sz = S — z is homeomorphic to IR^, and the concept of the winding number is well defined for Sz. 
Specifically, z can be considered as the point of infinity in M^. The winding number of x ^ w.r.t. C 
and z, denoted by wn(x; z, C) (C omitted when its choice is clear), is simply the summation of signed 
crossing numbers for any path connecting x to z. As in the planar case, we say that C is consistent w.r.t. z 
if wn(x; z, C) is either non-negative, or non-positive for all x G Sz. Similar to before, we define the total 
winding number w.r.t. a base point z as Tw(C; z) = Jg wn(x; z, C)dx. Let cr(P, Q; 0) denote the best 
cost to morph P to Q within domain Q.. 

Observation 5.2 If there is an optimal homotopy between P and Q that does not sweep through some point 
z G S, then we have cr(P, Q; S) = Q; Sz). 

Observation 5.3 Suppose H* is an optimal homotopy between P and Q with no anchor points. For any 
cell R in Arr(P + Q), ifH* sweeps through one point in its interior, then it sweeps through all points in R. 
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The simple proof for the above observation is in Appendix IE. II The key result is the following lemma, 
the proof of which can be found in Appendix IE.2I 

Lemma 5.4 If there is an optimal homotopy H* of P and Q with no anchor point, then the image of this 
optimal homotopy cannot cover all points in S. 

Given two homotopic paths P and Q from S sharing common endpoints, Lemma [54] and Obsei^vation 
I5.2l imply that if P can be morphed to Q optimally without anchor points, then there exists some point z G S 
such that (t(P, Q; S) = cr(P, Q; Sz). For this choice of z, it is necessary that the closed curve P o Q has 
consistent winding numbers in Sz. Once this z is identified, a{P, Q; Sz) is simply the total winding number 
of P o Q w.r.t. z, as suggested by Lemma 1431 because Sz is homeomorphic to the plane. Furthermore, by 
Observation 15.31 we only need to pick one point from each cell of Avv{P + Q) to check for the potential 
z. Specifically, let {zi, . . . , zi} be a set of such representatives, where / = 0{I). The optimal homotopy 
area (j{P, Q) is simply the smallest of all Tw(P o rev((5); Zj) for those ZjS with respect to whom the curve 
P o rev((5) has consistent winding numbers. Hence if we assume that if there is an optimal homotopy 
between P and Q with no anchor points, then we can compute cr(P, Q). 

Overview of the algorithm for sphere case. To compute the optimal homotopy between P and Q, 
we follow the same dynamic programming framework as before. If there is no anchor point in an optimal 
homotopy, then we use the discussion above to compute the optimal homotopy ai^ea. Otherwise, we identify 
the intersection point that serves as next anchor point, and recurse. The main difference lies in the component 
of computing <T(i,j) := o"(P[xj, Xj], (5[xj, Xj]), assuming that there is an optimal homotopy from P' = 
P[xj, Xj] to Q' = (5[xj, Xj] with no anchor points. Previously, this is done by checking whether P' o Q' has 
consistent winding numbers. Now, we need to check the same condition but against / = 0{I) number of 
potential representatives {zi, . . . , z;} as the potential point of infinity. This gives a linear-factor blow-up in 
the time complexity compared to the algorithm for the planar case. However, we show that this linear blow- 
up can be tamed down and we can again compute all a{r,j)s for all rs and all j > r in 0{P log/) time, 
after 0(n log n + A^) pre-processing time. See Appendix IE.3l for details. Overall, the total time complexity 
remains the same as before. 

Putting both cases (5 > and 5 = 0) together, we conclude with the following main result. 

Theorem 5.5 Given a triangulation K of an orientable compact 2-manifold M with genus g, let N be the 
complexity of K. Given two homotopic paths P and Q of total complexity n with I intersection points, we 
can compute an optimal homotopy and its area cr(P, Q; M) in 0{I^ log I + ng log n + A^) time. 

6 Conclusion 

In this paper, we propose a new cui-ve similarity measure which captures how hard it is to deform from one 
curve to the other based on the amount of total area swept. It is robust to noise (as it is area-based), and 
can be computed efficiently; to our knowledge, there is no other efficiently computable similarity measure 
for curves on surfaces. Our algorithm can be extended for cycles in the plane (see Appendix |Dll. It appears 
that our algorithm can also be extended to cycles on the surfaces. Indeed, if the optimal free homotopy 
has an anchor point, then we can break cycles into cui^ves that share a common start and end point, which 
then reduces to the problem of comparing curves on surfaces. However, on a surface the analog of Lemma 
ID. II no longer holds, so that two curves may intersect in M but not have an anchor point in the optimal 
homotopy; in this case, it is not immediately clear how to bound the size of the universal cover necessary for 
our algorithm. We leave the problem of working out these details, as well as improving the time complexity 
for comparing cycles, as an immediate future direction. 



15 



Currently, we assume that two input paths are simple paths which share staiting and ending points, which 
makes it easier to define homotopy equivalence. This leads to two natural questions, namely how to handle 
curves which do not share endpoints and how to deal with non-simple curves. Another interesting problem 
is to compute optimal isotopy area where we require that any intermediate curve during the deformation is 
also simple. 

Measuring similarity of curves on surfaces is an interesting problem, and many open areas remain. 
Geodesic Frechet-based measures ignore the topological constraints of underlying surface, while the ho- 
motopic Frechet distance, homotopy height, and our method require identification of a homotopy which 
optimizes some cost. As far as other measures of similarity which may be tractable, one interesting new 
idea would be to develop an area-based curve similarity measure that allows topological changes, such as 
allowing a region to be swept as long as it has trivial homology. Other directions include developing efficient 
curve simplification algorithms based on this measure, and studying similarity between curves from more 
general simplicial complexes than considered in this paper (such as a manifold with boundary or holes, or 
non-manifolds). 
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A Proof for Observation 3.1 



Note that H* is a map from □ M, where □ = [0, 1] x [0, 1] is the unit squai-e q i Q 

and a point (s, t) G □ will be mapped to (s). See the right figure for an illustration. 
(Since P and Q share starting and ending endpoint, the left and right sides of □ should 
be contracted to a point. We use the square view for simpler illustration.) The top and 
bottom boundary edges of this square are mapped to Q and P, respectively. Given an 
anchor point bj, let pj and be the parameters of bj in Hq and HI, respectively; that 
is, i^Q (pi) = ff^(q.j) = bj. By definition of anchor points, the pre-image of bj under — s ^' P 
the map H* necessarily includes a curve in □ connecting pi on the bottom edge to qj on the top boundary 
edge of □. Since bj / bj, the pre-images of bj cannot intersect with that of bj. Hence no two such curves 
can intersect each other, which means that pjS must be ordered in the same way as qjS. 



B Proof for Lemma 3.2 



Consider a time t in the homotopy, and let j = Ht. We first show that H deforms 7 consistently, so that 
every point on 7 is either fixed or deforms to the same side of 7. 

First note that if some portion of 7 is left sense-preserving at time t and then reverses its direction and 
becomes right sense preserving at time t+ a small amount later, some portion of the domain has been swept 
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(a) (b) (c) (d) (e) 



Figure 7: (a) and (b) p is fixed from 7 to 7+, but not so in 7++. Sweeping 7 to 7++ directly through the shaded 
region in (c) has a smaller area than first to 7+ then to 7++ (see shaded region in (d)). The darker shaded region in (d) 
is swept twice, (e) If the deformation changes orientation at 7, then there is a local fold in the regions swept. 




(a) (b) (c) (d) (e) 

Figure 8: (a) The cell R with highest positive winding number It boundary consists of alternating P-arcs (red) and 
Q-arcs (green). The two cases of relations between P[p, q] and R are shown in (b) and (d), respectively. For case (b), 
we can deform P to sweep through as shown in (c), and reduce the number of intersections by 2. Similarly, for case 
(d), we can identify any bigon R' and deform P to reduce the number of intersections by 2 as well. 

twice. Hence this homotopy cannot have minimal area, since we can create a smaller one which stops at 
time t and moves directly to some intermediate curve at a time greater than t+ without sweeping any portion 
twice. See Figure |7](e). 

Now suppose that some portion of 7 is deforming to one direction and another is morphing in the 
opposite direction. Since H is a homotopy and is therefore continuous, this means that there is at least one 
interval of fixed points between these two regions (which may possibly consist of a single point). Let p be 
an extremal point on this interval; see Figure 17] (a) for a picture when p is the only fixed point. In addition, 
since p is a fixed point but not an anchor point, where know there is some = t + dt where p is still on 

= 7+ and another t++ = t+ + (it where p is not on h^++ = 7++. 

Now we have several cases to consider. First, consider if H has directly reversed the direction of either 
portion of the curve (before p or after p), we are in a similar situation to the one previously discussed, since 
the curve goes from locally forward to locally backward (or vice versa). In this case, we again know that 
some area of the domain has been swept twice, which means 7++ has been swept over once and then was 
returned to, so we can reduce the area swept by H by reparameterizing H to move directly to 7++ without 
passing it and then reversing. (See Figure |7](b),(c), and (d) for an illustration.) 

Now if neither portion directly reverses, then 7++ must also be deforming to two different directions. 
Also, we know that 7++ must intersect 7 at some point q p, since we are essentially rotating around a 
central set of fixed points on these curves. In this case, we can again alter H to attain a smaller area swept 
by simply sweeping directly from 7 to 7''""*" ; this will reduce the area since the triangular region in the center 
bounded by 7, 7+, and 7++ will be swept only one time instead of twice. 

The claim thus follows, since any homotopy with no anchor points that is not sense preserving cannot 
be minimal. 
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C Proof for Lemma 



4.3 



We prove the claim by induction on the number of intersections between P and Q. The base case is when 
there is no intersection between P and Q. In this case, 7 is a Jordan curve which decomposes into two 
regions, one inside 7 and one unbound. By orienting 7 appropriately, every point in the bounded cell has 
winding number 1 and the claim follows. 

Now assume that the claim holds for cases with at most k — 1 intersections. We now prove it for the case 
with k intersections. Let an X-arc denote a subcurve of curve X. Consider the arrangement Arr(7) formed 
by 7 = Porev((5). Since P and Q are simple, every cell in this amngement has boundary edges alternating 
between P-arcs and Q-arcs. Assume without loss of generality that 7 has all non-negative winding numbers. 
Consider a cell R € Arr(7) with largest (and thus positive) winding number. Since its winding number is 
greater than all its neighbors, it is necessary that all boundary arcs are oriented consistently as shown in 
Figure [8](a), where the cell R (shaded region) lies to the right of its boundary arcs. 



If R has only two boundary arcs, e from P and e' from Q, respectively, then we 
can morph P to another simple curve P' by defomiing e through i? to —e' (where 
'— ' means reversing the orientation). See the right figure for an illustration. The p 
area swept by this deformation is exactly the area of cell R. Furthermore, after the 




Q 



deformation, every point x G R decreases their winding number by 1, and no other 



point changes its winding number. Since points in this cell initially has strictly positive 
winding number, the resulting curve j' = P' o Q still has all non-negative winding 




number. The number of intersections between P' and Q is k — 2. By induction hypothesis, cr(P', Q) = 
Tw(7'). Since Tw(7) - Tw(7') = Area(R), we have that Tw(7) = a{P', Q) + Area(R). It then follows 
from Observation ED and the fact a{P, Q) < a{P', Q) + Area(R) that a{P, Q) = Tw(7). 

Otherwise, the cell R has more than one P-arc. Take the P-arc ei with the smallest index along P, 
and let p be the ending endpoint of it. Let 62 be the next P-arc along the boundary of R, and q its starting 
endpoint, and Q[p, q] the Q-aic between ei and 62, denoted by e in Figured Obviously, the subcurve P[p, q] 
cannot intersect R, and P[p, q] and —Q\p, q] bound a simple polygon, which we denote by Q,. Either $7 is 
on the opposite side of the Q-arc e from the interior of R (Figure [8] (b)), or they ai^e on the same side (Figure 

EKd)). 

Case (1): R and Q are on the opposite side of e. In this case, the region is to the right of the oriented 
arc P\p,q]. Note that P does not intersect the interior of 0; as othewise, P will either intersect itself or 
intersect e, neither of which is possible. Hence only Q can intersect 0. Since Q is also a simple curve, there 
is no vertices of Arr(7) contained in the interior of il. As a result, every cell of Arr(7) contained in must 
have at least one boundary edge coming from P\p, q]. This implies that each cell contained in Q. has strictly 
positive winding number; that is, wn(a;; 7) > for any x € 17. This is because if a cell ^ C O has winding 
number 0, then its neighbor across its boundary on the other side of P[p, q] will have winding number — 1, 
as Q. is to the right of P[p, q]. This violates the condition that 7 has all non-negative winding numbers and 
thus cannot happen. 

We now deform P to P' by sweeping P[p,q] through Q. to Q\p,q]. See Figure [8] (c). The cost of 
this sweeping is Area(O) and Tw(7) — Tw(P' o Q) = Area(r2). P' is still simple, and the number 
of intersection points between P' and Q is now k — 2. Since wn(x; 7) > for any x G 17, we have 
wn(x; P' o Q) > for a; G $7. No other point will change their winding number after this deformation. 
Thus the curve P'oQ has all non-negative winding numbers as well. Hence by induction hypothesis, we have 
thata(P',Q) = Tw(P'oQ). Since o-(P, Q) -cj(P', Q) < Area(17) and Tw(7) -Tw(P'oQ) = Area(J7), 
it then follows from Observation I4.2l that o"(P, Q) = Tw(7). 

Case (2): R and $7 are both from the same side of e. We now consider the remaining case as shown 
in Figure [8] (d). Take the unbounded region 17 := IR^ \ 17 which is the complement of 17. This unbounded 
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region lies to the right of the oriented curve P[p, q\. Since both P and Q are simple, only Q can intersect 
$7 and there is no vertices of Arr(7) contained in the interior of Q. First, observe that it is not possible that 
$7 n (5 = 0. This is because otherwise, Q is the unbounded face of Arr(7) and thus the winding number 
for all points in is 0. This however is not possible as this will imply that any point y to the immediate 
left of P[p,q] has winding number —1, violating our assumption that all cells in Arr(7) have consistent 
(non-negative) winding numbers. 

Hence 17 n Q 7^ 0, and there are set of arcs from Q intersecting P\p, q]. Then there must exist a bigon 
cell R' bounded by only two arcs, one P-arc P[a, b] C P[p, q] and a Q-arc /3. See Figure [8] (d). Similar 
to the argument from the previous paragraph, we can show that points in R' must have strictly positive 
winding number. Now deform P to P' by sweeping P[a, b] through R' to /? as shown in Figure[8](e). P' 
is still simple, the number of intersection points between P' and Q is now k — 2. Only points in R' reduce 
their winding number by 1, and the resulting arrangement still has consistent winding numbers. As such, by 
induction hypothesis, we have that a{P', Q) = Tw(P' o Q). Since cr(P, Q) - a{P', Q) < Area(R') and 
Tw(7) - Tw(P' o Q) = Area(R'), it then follows from Observation SUthat a{P, Q) = Tw(7). 

D Cycles in the Plane 

We now consider the case where we have two simple cycles P and Q in the plane. We have the following 
characterization: 

Lemma D,l If the two simple cycles P and Q intersect, then there is an anchor point in the optimal homo- 
topy between them. 

Proof: Suppose that P and Q intersect but there is no anchor point in the optimal homotopy H* . By 
Lemma [3^ we know that H* must be sense preserving. However, this means that H* continually moves 
one curve to the other in one local direction, which means that one curve must be entirely contained within 
the other, contradicting the assumption that they intersect. ■ 

At this point, the algorithm for cycles which intersect each other reduces to the one for curves: if we 
know which intersection point between P and Q is the anchor point, we can simply "break" the cycles into 
two curves at this point; this will become the start and end point for each of the two curves. Hence our 
algorithm for cycles will take a multiplicative factor of 0{I) extra time than the algorithm for cui"ves, since 
we need to try each possible intersection point as the required anchor point. 

It then follows from Theorem I4.4l that: 

Corollary D.2 Given two polygonal cycles P and Q in the plane of n total complexity and with I > 
intersection points, we can compute the optimal homotopy and its area in 0{I{P log I + n log n)) time. 

The remaining case is that when the two polygonal cycles P and Q are disjoint. If one of the cycle 
contains the other, then the area of the optimal homotopy is simply the area sandwiched between these two 
simple cycles. This can be computed in 0{n log n) time easily. 

However, if the cycles aie disjoint but neither contains the other, then in a sense the "optimal" free 
homotopy between them will simply be the area bounded by each curve, since the homotopy can collapse 
each curve separately to a point and then deform the points to each other. Indeed, that the sum of area 
bounded by the two Jordan cycles P and Q is the minimum possible homotopy area follows from a similar 
argument as the proof of Observation 14.21 However, in this case, the free homotopy described above is not 
regular since it collapses a curve to a single point. Nevertheless, one can argue that there exists a sequence 
of regular homotopies whose areas converge to this sum. In other words, the optimal area homotopy is still 
well-defined (as the infinum of the area of regular- free homotopies between P and Q), although there does 
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not exist a regular homotopy to achieve this optimal area. (This is analogous to similar issues that arise in the 
general case, and is the reason for introducing more restricted integrals in the mathematical hterature ll28l .) 

E Missing Details for the Sphere Case 

Observation E.l Given a closed curve T and any two points z,w G S, we have that: wn(x;w) = 
wn(x; z) + wn(z; w) (all winding numbers are w.rt the curve T). In particular, for any two points zi, Z2 
from the same cell o/Arr(r), we have that wn(x; zi) = wn(a;; Z2) for all x 7^ zi, Z2. 

Proof: Let 7(2;, y) be a path connecting point x to y. Note that the concatenation between 7(x, z) and 
7(z, w) is a path from x to w. Since wn(x; w) is simply the summed signed crossing number of any path 
from x to w with respect to F, the claim follows immediately. ■ 

E. 1 Proof for Observations 15.31 

Proof of Observation 15.31 Suppose x and y are two points from the interior of R such that H* sweeps 
through x, but not y. Connect x with y by any path 7 in the interior of R. This path has to intersect the 
boundary of the region swept by H*, and let z be one such intersection point on 7. Obviously, there is a 
local fold in the optimal homotopy as it sweeps through z; namely, some intermediate curve will touch z and 
immediately trace back. Thus the input homotopy H* cannot be sense-preserving. Contradiction. Hence 
H* sweeps y as well. 

E.2 Proof for Lemma 

We prove the lemma by induction on the number of intersections between P and Q. When there is no 
intersection between P and Q (other than the common endpoints), the Jordan curve P o rev(Q) divides the 
sphere into two connected components, and the optimal homotopy is the smaller area of the two. The lemma 
holds for this base case. 

Now assume that the lemma holds for P and Q with at most k intersection points. We wish to show the 
result for the case where P and Q have k + 1 intersection points. Since H* has no anchor points, this optimal 
homotopy is sense-presei-ving by Lemma IJ!2] Assign an orientation to the closed curve C = P o Q so that 
locally, every point on the curve P will continuously deform to its right during the optimal homotopy. Now 
pick an arbitrary point z not on P and Q, and compute the winding number for each cell of Arr(P + Q) 
w.r.t. z. Take the cell R with the largest winding number. We assume that z ^ R. Suppose this is not the 
case and that z G i?. Then we show that we can change the choice of z to make this hold. 

Specifically, if z G ii, then the cell R must have winding number 0. Now take the cell R' of Arr(P + Q) 
with the smallest winding number, and let w be a point from R'. Obviously, wn(vif; z) < wn(x; z) < for 
any x G S. Now we consider the winding numbers w.r.t. to w instead of z. By Observation IE. 1 1 we have 
that wii(x; w) = wn(x; z) + wn(z; w). On the other hand, we have that wn(z; w) = — wn(w; z). Hence 
wn(z; w) > wn(x; z) > for any x G S. In other words, for this new choice of point w, we have that R 
still has the largest winding number and in this case, w ^ R. 

Similar to the proof of Lemma 14. 3 1 the boundary of this cell consists of alternating arcs from P and from 
Q, and they necessarily have the orientation as shown in Figure[T](a) (otherwise, one of the neighboring cell 
if R will have a larger winding number). Choose the P-arc ei that appears earliest along P, with p being its 
ending endpoint. Let P[q\ be the next intersection between P and R. We have that P[p, q] and e = Q[p, q] 
do not intersect each other. The Jordan curve P[p,q] o {rev{Q[p, q])) bounds two regions on the sphere 
(instead of a bounded one and an unbounded one in the planar case as shown in Figure [T](b) and (d)). We 
consider the region that lies to the right of P[p, q] (thus left of e), and denote it by Q. See Figure |9](a). Since 
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(a) 



(b) 



Figure 9: (a) Take Q as the region to the left of Q\p,q\. (b) There always exists a bigon 0' (possibly Q, is no 
Q-arc intersects P\p, q]. 

R is to the right of e, n = 0. As P is simple, only Q can potentially intersect the cell Q. Hence there 
always exists a bigon Q,' in which is a cell in Arr(P + Q). See Figure|9](b). Note that it is possible that 
il' = il. Let P[a, b] denote the P-arc that bounds the bigon Q,'. Let P' be a new curve obtained by replacing 
P[a, b] with (slightly above) Q[a, b]. Since P defomis always to the right in the optimal homotopy, and in 
the end, P[a,b] needs to defonn to some portion of Q (which is not necessarily Q[a,b] though), one can 
show that there is an optimal homotopy between P and Q that consists of first sweeping P[a, b] to Q[a, b] 
through Q,', and then optimally morph P' to Q. On the other hand, by the induction hypothesis, there is an 
optimal homotopy H' from P' to Q that does not sweep some point, say zi in S. There are now two cases: 

(i) If zi G S — il', then there is an optimal homotopy from P to Q that does not sweep zi as well. The 
induction step then holds and the claim follows. 

(ii) Otherwise, zi G 0'. Consider the cell P' G Arr(P' + Q) that contains zi. Note that P'n(S-O') / 0, 
as there is no vertices of Arr(P' + Q) contained neither on nor inside Q'. Hence R' must also contain 
some point, say Z2, that is outside of il'. It then follows from Observation 15 .3 1 that Z2 is not swept 
either. This leads us back to case (i), and the induction step again holds. 

The claim then follows by induction. 

E.3 Details of Algorithm for Sphere Case 

Dynamic programming framework. Similar to the planar case, let xq, . . . , x/ denote the intersection 
points between P and Q, ordered by their indices along P, with xq and x/ being the beginning and ending 
points of P and Q. Let T{i) denote the optimal homotopy area between P[0, x^] and Q[0, x^], and C[i, j] 
the closed curve formed by P[xj, Xj] o Q[xj, x^]. However, now we say that a pair of indices is valid 
as long as Xj and Xj have the same order along P and along Q. This is different from the definition of valid 
pairs of indices as in the planar case. 

Specifically, for any closed curve 7, it turns out that 7 can always have consistent winding number for 
some choices of the point of infinity z: that is, there always exists z G S such that wn(x, z; 7) is con- 
sistent for all X G S^. We call such choices of z consistent representatives w.r.t. 7. Let Tw*(7) := 
minz |Tw(7;Sz;)| where z ranges over all possible choices of consistent representatives w.r.t. 7. Then, 
Lemma |43] and [53] imply that if there is an optimal homotopy between P[xi,Xj] and Q[xj,Xj] with 
no anchor points, then the optimal homotopy area is Tw*(7). However, different from the planar case, 
Tw*(C[i, j]) is defined for all valid pairs of i, js, and it may not in general be the optimal homotopy area 
for P[xi, Xj] and (5[xj, Xj]. We now have the following recurrence: 




and (j, i) is vahd i Tw*(C[i, i]) + T{j) }, otherwise 



if i == 
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As before, the final goal is to compute T{I) = a{P, Q). 



Q / , *J 






** 







Computing Tw*s. Here we describe how to compute Tw*(C[i, j]) efficiently. Specifically, we show how 
to compute all Tw*(C[r, j])s for all j > r in 0{I) time, for any fixed r, after 0{n\ogn + /log/ + N) = 
0{n log n + A^) preprocessing. 

First, let us choose the representatives by taking two points around each 
intersection points Xj between P and Q as shown in the right figure. Consider 
only those representatives to the right of Q (which are those above Q in the 
right figure), and denote them by Z = {zi , . . . , z/}. Z is sorted by their indices 
along Q. Those to the left of it can be handled similarly. The first obsei^vation 
is that for any two consecutive representatives, wn(3;;zj) — wn(x;zj+i) is 1 or —1, depending on the 
orientation of P-arc that separating them. 

Now to compute which Z j will give consistent winding numbers, we first compute the winding number 
of each cell in Arr(/' + Q) for zi. Next, take the cells Ri and R2 with minimum and maximum winding 
numbers, and assume that Zjj and Zjj are their representatives. If there are more than one cells with largest 
(or smallest) winding numbers, just pick an arbitrary one. By Observation IE. 1 [ the closed curve C[r, j] := 
P[xr,Xj] o (5[xr,Xj] has all non-negative winding number w.r.t Zj^, and all non-positive winding number 
w.r.t. Zjj. Hence we simply compute the total winding number Tw(C[r, j]; Zj^ ) and Tw(C[r, j]; Zjj), and 
return the one with the smaller absolute value as Tw* (C[r,j]). We refer to the indices ii and 12 as the wn- 
min and wn-max indices, respectively, and these two total winding numbers as valid total winding numbers. 
Basically, by Lemma lE!2l below. the smaller of the absolute values of the two valid total winding numbers is 
simply the best cost to deform P[xr, Xj] to Q[Xf., Xj] without using anchor points. This improves the time 
complexity of computing each Tw*(C[r, j]) to 0(/) time, instead of the naive 0{ln) time by computing 
all Tw(C[r, j]; Zj)s, for i G [1, /], from scratch. 

Lemma E.2 Given an arbitrary oriented (not necessarily simple) curve C = P' o rev{Q') on S, let 
Arr(C) be the arrangement of C, and Z = {zi, . . . ,2;^.} a set of representative points from each cell in 
C. Pick an arbitrary point, say z\, and compute the winding number of C w.r.t. each Zi. Let ii and 
12 be the wn-min and wn-max indices. Then the best cost to deform P' to Q' with no anchor point is 
mm{\Tw{C;zi,)\,\Tw{C;z,,)\}. 

Proof: First, call a point z valid if wn(x; z, C) is consistent for all x E S. The optimal cost to deform P' 
to Q' with no anchor point is minyaiid zez |Tw(C; z)\. For any base point z, note that by Observation IE. II 
we have wn(x; z) = wn(a;; zi) — wn(z; zi). In order for the winding number to be consistent, we need that 
either wn(x; 2;) > for any x G S, or wn(x; 2;) < for any x G S. Assume it is the former case. Then 
wn(2;; zi) < wn(x; zi) for all x G S, implying that wn(2;; zi) = wn{zi-^ ; zi). Furthermore, note that 



Tw(C;2;)= / wn{x; z)dx = / [wn{x; zi)—wn{z; zi)]dx = / [wn(x; zi)— wn(2;jj ; zi)]dx = Tw(C; Zj^) 
Js Js Js 

If it is the latter case, then wn(z; zi) > wn(x; zi) for all x G S, implying that wn(z; zi) = wn(zj2 ; zi). In 
this case we have that 

Tw(C;z)= / wn{x ; z)dx = / [wn(x; zi)— wn(z; zi)](ix = / [wn(x; zi)— wn(zj2; zi)]dx = Tw(C; z^j) 
Js Js Js 

The optimal cost minyaiidzez |Tw(C;z)| is thus achieved by the smaller one of the absolute value of 
Tw(C; Zi,) and Tw(C; ZiJ. ■ 

To further improve the time complexity, we will start with C[r, r + 1], and update the winding number 
information in each cell as well as the valid total winding numbers, as we traverse P and pass through each 
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intersection point Xj. To this end, we use the same range tree data structure as in Section l42l Specifically, 
we use this data structure to maintain the winding number information w.r.t a fixed based point zi. The 
wn-min and wn-max indices can be easily maintained by storing at each internal node the minimum and 
maximum winding number within its subtree. We can also maintain the total winding number w.r.t. the base 
point zi- The time complexity for updates remains the same as before (i.e, 0(log /) time per update). 

The remaining question is to compute the valid total winding numbers as i increases. Let A denote the 
total area of topological sphere S. First, observe that for a fixec curve C, by Observation lE.ll we have 
wn(3;; z, C) = wn(x; zi, C) — wn(z; zi, C) with respect the fixed based point zi. Hence we have that: 



Assume that ii is the wn-min index and 12 is the wn-max index. Hence we can compute Tw(C;zjJ 
and Tw(C;zj2) in 0(1) time using Eqn ([T]), since Tw(C;zi), wn(zj^;zi,C) and wn(zj2;z2,C) are all 
maintained as C changes from C[r, u] to C[r,u + 1]. 

Putting everything together, with 0(n log n + N) pre-processing time, we can compute all a{r,j)s for 
all j > r for any fixed r, in 0(/log/) time, and thus computing all Tw*(C[r, u])s for all r G [1,/] and 
all u < r,m 0{P log I) total time. Putting everything together, the dynamic programming problem can be 
solved in 0{n log n + l"^ log I + N) total time. 




Tw(C;zi) - ^•wn(z;zi,C). 



(1) 
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